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The known galaxies most dominated by dark matter (Draco, Ursa Minor and 
Andromeda IX) are satellites of the Milky Way and the Andromeda galaxies 1-4 . 
They are members of a class of faint galaxies, devoid of gas, known as dwarf 
spheroidals 3 5 , and have by far the highest ratio of dark to luminous matter 3 ' 6 . 
None of the models proposed to unravel their origin 7 10 can simultaneously 
explain their exceptional dark matter content and their proximity to a much 
larger galaxy. Here we report simulations showing that the progenitors of these 
galaxies were probably gas-dominated dwarf galaxies that became satellites 
of a larger galaxy earlier than the other dwarf spheroidals. We find that a 
combination of tidal shocks and ram pressure swept away the entire gas content 
of such progenitors about ten billion years ago because heating by the cosmic 
ultraviolet background kept the gas loosely bound: a tiny stellar component 
embedded in a relatively massive dark halo survived until today. All luminous 
galaxies should be surrounded by a few extremely dark-matter-dominated dwarf 
spheroidal satellites, and these should have the shortest orbital periods among 
dwarf spheroidals because they were accreted early. 



Draco, Ursa Minor and Andromeda IX have mass-to-light ratios (M/L) larger than 100, 
but the majority of the other dwarf spheroidals (dSphs) in the Local Group have a lower 
M/L, of order 10-30 1 , typical among dwarf galaxies 11 ' 12 . Another important difference 
is that Draco and Ursa Minor nearly stopped forming stars more than 10 billion years 
ago, while other dSphs continued to form stars for many billions of years 2 . The modest 
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potential well of these extreme dwarfs cannot be the one property that determined their 
nature. Their halo masses are too large to invoke suppression of gas accretion owing 
to the cosmic ultraviolet background at high redshift 9 ' 13 or blow-out due to supernovae 
winds 14 . Tidal shocks occurring as a dwarf repeatedly approaches the primary galaxy can 
transform rotationally supported systems resembling dwarf irregular galaxies (dlrrs) into 
systems dominated by random motions, similar to dSphs 10 . This tidal stirring can explain 
why dSphs are more clustered around the primary galaxies relative to dlrrs but it leaves a 
significant gas component inside the dwarf, so that star formation can continue for several 
billions of years instead of being truncated early. Ram pressure in a hot gaseous corona 
could strip their gas completely 7,15 but the limitations of simulations so far have not allowed 
for firm predictions. For instance, existing calculations keep the structure of the stars and 
halo fixed in time and neglect radiative cooling and heating of the gas 16,17 . 

These earlier studies have explored the effect of a single gas removal mechanism and 
are, at best, only qualitatively consistent with the current structure formation paradigm, 
a model with cold dark matter and a cosmological constant (ACDM). Recent attempts 
to study the evolution of dwarf satellites directly in cosmological simulations rely on 
semi-analytical methods to model the baryonic component rather than solving the fluid 
equations 18 . These models neglect ram pressure, and since stripping by tides is slow and 
inefficient, they cannot explain the complete absence of gas and early truncation of star 
formation of the darkest dSphs. 

In cold dark matter models the present-day spatial distribution of subhalos within 
primary halos retains some memory of their infall time 19 . Satellites orbiting closer to 
the primaries were on average accreted earlier than those orbiting at larger distances. 
Interestingly, Draco and Ursa Minor lie at 68 and 86 kpc, respectively from the Milky Way, 
Andromeda IX at 45 kpc from Andromeda, while other dSphs orbit as far as 200 kpc from 
the primaries 1 . We use a high resolution ACDM dark matter-only cosmological simulation 
of the formation of a Milky Way-sized halo 20 (see also Supplementary Information). At 
z = we identify three subhalos having distances below 100 kpc from the center and with 
a peak circular velocity in the range 25 — 30 km/s 6 (Figure 1). We track the orbits of the 
satellites back in time and find that two of them were accreted early, between z = 2.5 and 
z = 1.5. 

We then construct a high resolution N-body + smoothed particle hydrodynamics (SPH) 
model of a dwarf galaxy satellite having a disk of stars and gas inside a cold dark matter halo 
(Figure 1) with a peak velocity of about 40 km/s, comparable to that of the two identified 
cosmological subhalos before they were accreted onto the Milky Way (see Supplementary 
Information). We assume that 80% of the baryonic disk mass is in a gas component. The 
inefficient conversion of gas into stars is expected at these low mass scales since most of 
the gas will have densities below the threshold for star formation 21 . In addition, at z > 2 
the gas in the dwarf is heated to a temperature of over 10 4 K and ionized by the cosmic 
UV radiation, which further suppresses star formation. 

The dwarf model is placed on an eccentric orbit inside a massive Milky Way-sized halo 
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Figure 1. Morphological evolution of the dwarf satellite galaxy. Left: Color-coded 
projected density map of the cosmological run at z = 0; the box is 260 kpc on a side, which 
corresponds to the virial radius of the Milky Way-sized halo. The peak density along the line 
of sight is shown, ranging from 10~ 29 — 10~ 24 g/cm 3 , with the color coding from blue (lowest 
density), through yellow, then red, to brown (highest density). The three satellites that meet the 
distance and circular velocity constraints (see text) are highlighted with black boxes. Right: From 
top to bottom the insets show the stellar component of the dwarf galaxy, color coded in projected 
density, at different times. Only regions with densities in the range 10 -28 — 10~ 23 g/cm 3 are 
shown, with the color coding from dark red (lowest density) to yellow (highest densities). In the 
top inset, the initial disk is shown edge-on (the box is 8 kpc on a side). In the left middle inset, 
the system is close to the second pericenter passage; the stars have assumed a strong bar-like 
configuration and heating is evident in the outskirts (the box is 7 kpc on a side). In the bottom 
inset, the end state is shown; the bar has been heated into a diffuse spheroid and any disk-like 
feature has been erased (the box is 4 kpc on a side and a a projection along a random line of sight 
is shown). In the right middle inset, the trail of gas produced by ram pressure is shown, while 
even the residual gas in the center is stripped. The color-coded gas density projected onto the 
orbital plane is shown (densities in the range 1CP 30 — 1CP 23 g/cm 3 , color coding as above), and 
the box is 100 kpc on a side. 
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model which is a replica of that in the cosmological simulation. We include radiative 
cooling as well as the heating and ionizing flux from the cosmic ultraviolet background 
radiation 22 , and we embed a diffuse gaseous halo inside the dark halo of the primary. Such 
a halo is expected as a by-product of the process of galaxy formation and has a density and 
temperature consistent with observational constraints (see Supplementary Information). 

With an orbital time of about 1.7 Gyr the dwarf undergoes as many as 5 pericenter 
passages in 10 Gyr. At the first pericenter passage its dark halo loses 60% of its mass. 
The disk, deep inside the potential well of the halo, suffers no stripping, but the tidal 
perturbation triggers a strong bar instability (Figure 1) and simultaneously heats the stars 
in the disk. The bar funnels most of the gas towards the central kiloparsec. Gas will be 
removed from the galaxy if the ram pressure force exerted by the diffuse hot gas, which 
is proportional to p g V 2 (p g being the density of the gaseous halo and V being the orbital 
speed of the galaxy), exceeds the gravitational restoring force provided by the potential 
well of the galaxy 23 . Ram pressure readily removes the gas outside the bar radius on 
a timescale of less than 10 8 years, but not the more tightly bound gas inside the bar. 
When the satellite crosses the pericenter a second time, the new tidal shock lowers the 
halo density by a factor of 2 inside 1 kiloparsec, so that Vpeak drops to less than 30 km/s. 
Since the potential well has become shallower, even the gas sitting inside the bar can be 
swept away by ram pressure (Figure 1, see also Supplementary Information). The cosmic 
UV background heats and ionizes the gas, which enhances stripping significantly since the 
higher gas pressure opposes the gravitational restoring force. Once the first two orbits have 
been completed no gas is retained by the dwarf. 

The response of the system to the tides becomes progressively more impulsive at each 
new pericenter passage. As a result, the initially disky stellar distribution is heated into a 
nearly spherical, isotropic configuration (Figure 1). The expansion reflects the attempt of 
the system to gain a new equilibrium as the internal binding energy is lowered. The removal 
of gas due to ram pressure is crucial for tidal heating to be effective (see Supplementary 
Information). 

After 10 billion years a diffuse spheroidal galaxy has replaced the gas-rich disk. During 
the last few orbits the stellar velocity dispersion profile is fairly flat (left panel of Figure 2), 
matching the profiles observed for Draco and Ursa Minor 24 . The initial angular momen- 
tum of the stars has been transported outward during the morphological transformation 8 , 
producing a ratio between rotational velocity and random motions v rot /cr < 0.2 in the 
inner few hundred parsecs. The surface brightness and total luminosity resemble those of 
the faintest dSphs (see caption of Figure 2). A substantial halo is preserved within a few 
kiloparsecs from the center so that the final Mj L is larger than 100 (right panel of Figure 
2). The central dark matter density, which has decreased by almost a factor of 4 since 
the beginning of the simulation, is ~ O.2M /pc 3 , comparable to that of Draco and Ursa 
Minor 1 . 

We predict that all massive galaxies should have a few extremely dark matter dominated 
satellites as the mechanism reported here is completely general within hierarchical structure 
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Figure 2. Structural properties of the simulated dwarf after 10 billion years of evo- 
lution. Left: The line-of-sight stellar velocity dispersion profiles are shown for three random 
directions (black lines) perpendicular to each other, together with published observational data 
points for Draco (red triangles) and Ursa Minor (green squares), including formal l-cr error bars 
24 . The curves are shown out to the radius for which data points are available. Right: Mass 
profiles shown out to the radius (from the center of the simulated dwarf) at which stars are grav- 
itationally bound. The dashed lines show the initial (black dashed line) and final (red dashed 
line) ratios of the total mass to the baryonic mass. The solid lines show the initial (black solid 
line) and final (red solid line) B band mass-to-light ratios of the dwarf. We have assumed stellar 
mass-to- light ratio (M/Lb)* = 1.5 (initial) and (M/Lb)* = 5 (final) to compute the initial and 
final luminosities, Mb = —12.5 and Mb = —9. The final central surface brightness is [ib ~ 26 
mag arcsec -2 . We note that an initial (M/Lb)* = 1-5 is motivated by the fact that at z > 2 a 
stellar population is at most three billion years old, whereas a final (M/Lb)* = 5 is consistent 
with passive fading of the stellar population for about ten billion years 10 . 



formation. The efficiency of star formation in isolated low-mass galaxies can be lower than 
that assumed here 21 (see Supplementary Information), and hence dwarf satellites having 
gas fractions higher than 80% of the baryons before stripping probably did exist. Once 
they lose their gas, these systems will turn into dSphs even more dark matter dominated 
than Draco and Ursa Minor. They should orbit the Milky Way and M31 and have escaped 
detection so far owing to their low surface brightness, especially at low Galactic latitude 25 . 
Their halo masses should be comparable to those of known dSphs and thus might help to 
solve the missing satellites problem 26 . The Ursa Major dwarf recently detected by the Sloan 
Digital Sky survey 27 ' 28 ) which has a halo as massive as that of Draco but a luminosity a 
hundred times smaller, is probably one of them. Three more systems such as Ursa Major 
would be enough to bring theory and observations in agreement at the high end of the mass 
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function of satellites, corresponding to a peak circular velocity of around 25 — 30 km/s, 
where suppression of baryonic infall by reionization is hardly effective. 

When the ultraviolet background is an order of magnitude weaker, as is predicted 22 at 
z < 1, about 30% of the original gas component stays in the galaxy (see Supplementary 
Information). Therefore dwarfs that fell into the Milky Way halo late can continue forming 
stars, and tidal shocks will produce periodic bursts of star formation 10 . These newcomers 
should account for those dSphs that have fairly normal mass-to-light ratios, extended star 
formation histories and larger distances from the primaries 1 ' 2 . This explains why Fornax is 
ten times brighter than Draco and has a very different star formation history despite having 
a comparable depth of the potential well 1 ' 6 ' 29 . It implies that there should be a positive 
correlation between M/L and the infall epoch of dwarfs, and thus a negative correlation 
between M/L and their orbital time. The dSph Tucana represents the biggest challenge 
to our model because it lies far from any massive galaxy (see Supplementary Information). 
The accurate determination of the orbits of the dwarfs expected from on going and future 
astrometric missions such as the Space Inteferometry Mission and the Global Astrometric 
Interferometer for Astrophysics will be able to test this prediction. 
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SUPPLEMENTARY INFORMATION 

Here we briefly describe the setup of the initial conditions and the numerical methods used 
to perform the simulations presented in the Letter. This is followed by a discussion of the 
tests and the analysis performed in order to verify and understand the results in depth. 

1 Simulation code 

We have used the fully parallel, N-body+smoothed particle hydrodynamics (SPH) code 
GASOLINE to compute the evolution of both the collisionless and dissipative component 
in the simulations. A detailed description of the code is available in the literature 1 . Here 
we recall its essential features. GASOLINE computes gravitational forces using a tree- 
code 2 that employs multipole expansions to approximate the gravitational acceleration on 
each particle. A tree is built with each node storing its multipole moments. Each node is 
recursively divided into smaller subvolumes until the final leaf nodes are reached. Starting 
from the root node and moving level by level toward the leaves of the tree, we obtain a 
progressively more detailed representation of the underlying mass distribution. In calcu- 
lating the force on a particle, we can tolerate a cruder representation of the more distant 
particles leading to an O(NlogN) method. Since we only need a crude representation for 
distant mass, the concept of "computational locality" translates directly to spatial locality 
and leads to a natural domain decomposition. Time integration is carried out using the 
leapfrog method, which is a second-order symplectic integrator requiring only one costly 
force evaluation per timestep and only one copy of the physical state of the system. In 
cosmological simulations periodic boundary conditions are mandatory. GASOLINE uses a 
generalized Ewald method 3 to arbitrary order, implemented through hexadecapole. 

SPH is a technique of using particles to integrate fluid elements representing gas 4 ' 5 . 
GASOLINE is fully Lagrangian, spatially and temporally adaptive and efficient for large 
N. The version of the code used in this Letter includes radiative cooling and accounts 
for the effect of a uniform background radiation field on the ionization and excitation 
state of the gas. The cosmic ultraviolet background is implemented using the Haardt- 
Madau model 7 which includes photoionizing and photoheating rates produced by QSOs 
and galaxies starting at z — 7. The assumption that reionization occurred at z = 7 is 
conservative but consistent with the combination of the 3rd year WMAP results and the 
Gunn-Peterson effect in the spectra of distant quasars 8 . We use a standard cooling function 
for a primordial mixture of atomic hydrogen and helium (the metallicity in dSphs is indeed 
much lower than solar 9 , with — 1 < [Fe/H] < —2). The cooling shuts off below 10 4 K. The 
internal energy of the gas is integrated using the asymmetric formulation, that gives results 
very close to the entropy conserving formulation 6 but conserves energy better. Dissipation 
in shocks is modeled using the quadratic term of the standard Monaghan artificial viscosity 5 . 
The Balsara correction term is used to reduce unwanted shear viscosity 10 . 
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2 Initial Conditions: the cosmological simulation 

The cosmological simulation employs the so called volume renormalization technique to 
achieve higher resolution in a region of interest (described in the next section). It is one 
of the simulations already presented elsewhere 11,16 . The virial mass of the halo that was 
selected at z = to be re-simulated at higher resolution is 1.15 x 1O 12 M (the virial mass is 
measured within the virial radius R v ir, the radius enclosing an overdensity of 100 times the 
critical density p C rit)- A recent analysis of all observational constraints points to a Milky 
Way halo of about 10 12 M o n . The halo was originally selected within a low resolution, 
dark matter only simulation run in a concordance, flat, A-dominated cosmology: Q = 0.3, 
A=0.7, h = 0.7, <t 8 = 0.9, shape parameter T = 0.21, and tt b = 0.039 12 . The size of the 
box, 100 Mpc, is large enough to provide realistic torques. The power spectra to model the 
initial linear density field were calculated using the CMBFAST code to generate transfer 
functions. Dark matter particle masses in the high resolution regions are 8.05 x 10 5 M , 
and and the force resolution, i.e. the gravitational softening, is 0.3 kpc. In total there are 
1.4 x 10 6 particles within the virial radius. 

With our choices of particle number and softening, the smallest subhalos resolved have 
typical circular velocities of 10% of their host. For all particle species in the high resolution 
region, the gravitational spline softening, e(z), evolved in a comoving manner from the 
starting redshift (z ~ 100) until z=9, and remained fixed at its final value from z=9 to 
the present. The softening values chosen are a good compromise between reducing two 
body relaxation and ensuring that disk scale lengths and the central part of dark matter 
halos will be spatially resolved. Integration parameter values were chosen according to the 
results of previous systematic parameter studies 14 . 

The merging histories and angular momentum of the parent dark matter halo are sup- 
posed to play a major role in defining the final properties of the galaxies that form inside 
them. It is then important to make sure that our halos have merging histories and spin 
parameters somewhat representative of the global population. The halo was selected with 
the only criteria that the redshift of the last major merger (z/ mm ) was > 2 and that there 
are no halos of similar or larger mass within a few virial radii (a major merger is defined 
here as having a 3:1 mass ratio). The halo spin parameter, A, is ~ 0.05 at z — 0, close to 
the average value for cosmic halos ~ 0.035 15 . Cosmological simulations with hydrodynam- 
ics that use this same and other similar initial conditions show that realistic disk galaxies 
similar to the Milky Way do arise in these dark matter halos 16 . 

2.1 Volume renormalization technique 

The large dynamic range involved in cosmological simulations aimed at resolving the scale 
of galactic halos calls for techniques that concentrate the computational power on the 
object of interest. This is achieved by the volume renormalization technique 17 . A large 
scale simulation is done at low resolution and the regions of interest are identified, e.g. 
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QSOs forming at high redshift, a large cluster of galaxies or a Milky Way-sized halo such 
as in this Letter. Next, initial conditions are reconstructed using the same low-frequency 
waves present in the low resolution simulation but adding the higher spatial frequencies. 
To reduce the number of particles and make the nonlinear simulation possible with the 
same cosmo logical context, we construct another set of initial conditions with particles 
whose mass, and therefore mean separation, increase with the distance from the center of 
our volume. Note that because tides are important in the formation of the filaments it is 
not sufficient to extract just the central region. Generating initial conditions (ICs) in this 
approach is nontrivial. The ICs are calculated on a grid determined by the initial mass 
resolution. The power spectrum is realized by using Fast Fourier Transforms (FFTs) to 
determine displacements of particles from this grid. We use a particular technique that 
allows one to calculate high resolution FFTs only in the regions of the simulations were 
mass and force resolution need to be high. This method has been successfully used in a 
wide range of cosmological studies 18 ' 19 ' 20 ' 11 . 

3 Initial Conditions: the dwarf galaxy model 

The multi-component dwarf galaxy models are built using a standard technique that pro- 
duces initial conditions very close to equilibrium 20 ' 21 . The model is further evolved in 
isolation for several Gyr in order to remove possible transients caused by initial particle 
noise, and the end state is used as initial condition. The scheme we use to assign the 
structural parameters of the baryons in a dark matter halo assumes that baryons collapse 
into the halo having the same specific angular momentum distribution of the halo and con- 
serving angular momentum during collapse 23 . The properties of the baryonic component 
are related to those of the dark halo, namely its mass, angular momentum (via the spin 
parameter) and concentration. We refer to previous papers for a detailed description of the 
scaling relations that relate the baryonic component and the dark halo 24 ' 25 . The technique 
has been widely used to build models of dwarf galaxies 26 ' 27 . 

The bayonic component comprises an exponential disk of gas and stars embedded in 
an NFW halo 28 . The halo has Vp Ca k = 40 km/s. In CDM models low mass halos have on 
average higher concentration parameters c = Ryh/ r s (where i? vir is the virial radius and r s 
the scale radius) relative to higher mass halos and, for a given halo mass, the concentration 
decreases with redshift roughly as (1 + z)~ l . We choose c = 20 for the halo of the dwarf. 
This value of the halo concentration matches the average concentration of the subhalos 
with Vp Ca k in the range 30 — 40 km/s in the cosmological simulation before they fall into 
the primary halo and is consistent with the typical values expected at these mass scales 
at z ~ 2, as shown elsewhere 29 . In choosing the radius and mass of the disk there is 
some freedom. However, there are constraints coming from the link between disk radius 
and the spin parameter of the halo in which the disk forms, and from the information of 
the typical fraction of the mass of present-day dwarf galaxies that is contributed by the 
baryonic disk. The first constraint brackets the range of possible disk radii at given disk 
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Figure 3. Structural properties of the initial dwarf galaxy model. The left panel shows the 
initial rotation curve. The red solid line shows the contribution of the gas, the red dashed line 
the contribution of the stars, the black dashed line shows the dark matter contribution and the 
black solid line shows the total rotation curve resulting from the sum of the contribution of the 
different components. The right panel shows the initial baryonic surface density profile. 

mass, the second the range of possible disk masses for a given halo mass. Since our simple 
working hypothesis is that the progenitors of dwarf spheroidals had baryon fractions in the 
range of the normal values of galaxies in the present-day Universe, we fix the disk mass 
fraction to fd = 0.035M vir . This value is typical in mass models that reproduce dwarf and 
low surface brightness galaxies 29 and is consistent with the baryonic Tully-Fisher relation 
for galaxies with rotational speeds in the range 40 — 60 km/s 30 . The exponential disk scale 
length is then fixed by the adopted spin parameter, A = 0.05. The log-normal distribution 
of spin parameters of cosmological halos peaks at around 0.035 15 but (1) 40% of the halos 
have larger spins and (2) modeling of their rotation curves suggests that dwarf galaxies 
have an average spin larger than the mean value of the galaxy population as a whole 29 . 
The resulting exponential disk scale length, which is a function of A, the halo concentration 
c and /d 24 , is 0.55 kpc. 

The chosen disk mass and radius yield an initial central baryonic surface density ~ 
4M /pc 2 . We choose the stellar mass fraction in such a way that it is consistent with the 
amount of star formation expected in a galaxy with the assigned structural parameters. 
Within the standard structure formation paradigm galaxies in halos of 10 10 M Q or less, 
namely in the range of our model, should have little or no star formation because their 
disks are gravitationaliy stable nearly everywhere 31 , having Toomre parameters Q > 1. The 
stability arises because lower mass halos should harbor lighter disks. The latter statement 
is true if the disk mass fraction either varies little or tends to decrease for decreasing 
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halo masses, in agreement with the analysis of galaxy rotation curves across a range of 
masses 30 ' 31 . Observations show examples of galaxies that are Toomre-stable and yet have 
on-going or recent star formation, albeit at a fairly low level 32 ' 33 . Recent models can account 
even for such departures from the Toomre criterion 34 ' 35 . In one such model the threshold 
for star formation is determined by the minimum surface density required for molecular gas 
to form 34 . In our dwarf galaxy model about 1.5 x 10 7 M Q of baryons have surface densities 
above the predicted critical surface density predicted, which is ~ 2M /pc 2 . We thus assign 
to our model a stellar mass of 1.5 x 1O 7 M , corresponding to only 20% of the baryons. Since 
the critical surface density for star formation can rise by up to a factor of 5 in presence of 
ionizing radiation fields such as the cosmic ultraviolet background 34 , the assumed stellar 
mass fraction should be viewed as an upper limit for a reasonable progenitor at these 
mass scales. Although the criterion based on surface density is convenient, it is known to 
have difficulties in reproducing some observations. Therefore we compared our choice with 
another, more recent work that predicts the stellar fraction using fully self-consistent N- 
body+SPH simulations with a physically motivated star formation recipe and sink particles 
to represent sites where molecular clouds and star clusters would form 35 . These simulations 
employ models of galactic disks with a variety of initial structural parameters, gas mass 
fractions and gas temperatures. They explore the dependence of the star formation rate on 
the initial minimum Toomre parameter of the combined gaseous and stellar disk, Q sg , and 
find that star formation is possible up to values of Q sg ~ 1.6, significantly larger than the 
standard threshold Q g ~ 1, where Q g is the local Toomre parameter of the gas. Moreover, 
by running a large set of simulations for about 2 Gyr, they derive a relation between the 
minimum initial Toomre parameters of the systems, Q sg (min) or Q g (min), and the total 
amount of gas that can be converted into stars. We compared our initial setup with their 
results. Assuming a fully gaseous disk at the start, our initial galaxy model would have 
Q sg (min) = Q g (min) ~ 1.42 (reached between 1 and 2 disk scale lengths). For such a 
value of Q g (min), Figure 15 of their paper (bottom panel) indicates that between 10% and 
17% of the initial gas should be converted into stars (13% if we compare with their best 
fitting curve). These fractions are close to the 20% that we assumed, though somewhat 
smaller; indeed, as we emphasize in the Letter, a smaller initial stellar fraction would make 
our scenario even stronger. The rotation curve and initial baryonic surface density profile 
of the dwarf galaxy model are shown in Figure 3. 

We choose the mass and force resolution in the dwarf galaxy models in such a way that 
possible numerical artifacts are under control. We adopt 2 x 10 5 particles in the gaseous disk 
and 10 6 particles in the stellar disk and dark matter halo, for a total of 2.2 x 10 6 particles. 
The spatial resolution of the dark halo and stars is set by the gravitational softening which is 
equal to 35 pc. The spatial resolution in the gas component is given by the SPH smoothing 
length, defined as the radius of the spherical kernel adopted by GASOLINE, which contains 
32 particles. The smoothing length is adaptive in space and time. In the initial conditions 
it is about 25 pc at about one disk scale length, and it drops to ~ 15 pc during the stage 
of the bar driven gas inflow close to the first pericenter passage. Since most of the mass in 
an exponential disk is contained within about 3 scale lengths, or 1.5 kpc, the gaseous disk 
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is resolved by about 60 resolution elements initially; later, as the bar forms, the average 
smoothing length decreases, but also the size of the gas component becomes smaller by 
almost a factor of 2 (due to the gas inflow and stripping of the outer regions), so that 
effectively the resolution remains close to 60 smoothing kernels. This resolution is above 
that required to properly model ram pressure stripping based on previous works which 
found that 50 resolution elements in a linear dimension (or 25 along a radius) is needed 36 ' 37 . 
Moreover, both the smoothing length and the gravitational softening of the gas are at all 
times much smaller than the local Jeans length, and have comparable values. This ensures 
that gravity and pressure forces are correctly balanced at the Jeans length, thus avoiding 
spurious fragmentation or artificial suppression of collapse 38 . The Jeans length is defined 



as Aj = y§^) where v s is the sound speed, G is the gravitational constant and po is the 

(local) density. It is smallest during bar formation since the temperature of the bound 
gas is ~ 2 x 10 4 K at all times while the density peaks at this time, being close to 10~ 23 
g/cm 3 near the center. Yet even in this extreme situation Aj ~ 300 pc, i.e. about an 
order of magnitude larger than either twice the SPH smoothing length or the gravitational 
softening 38 . 

4 Initial Conditions: the Milky Way halo model 

The live dark halo model used for the Milky Way is based on the Milky- Way sized halo 
of the cosmological simulation (see above). We use a slightly larger value of the circular 
velocity at the virial radius, V vir = 180 km/s, in order to have Vp ea k = 230 km/s at about 
10 kpc from the center without including a baryonic disk component. The concentration 
is c = 10. With our choice of parameters the shape of the rotation curve is quite flat 
out to 100 kpc and the resulting halo mass within 100 kpc is ~ 8 x lO n M , consistent 
with the orbital dynamics of the Magellanic Clouds and distant dSphs 39 ' 40 . We use 10 6 
particles to model the halo and a softening of 300 pc, as in the cosmological simulation. 
We embed a gaseous halo in hydrostatic equilibrium within the live dark matter halo. 
The density of the gaseous halo near the pericenter of the orbit of the satellite, 30 kpc, is 
2 x 10~ 4 cm -3 , comparable to that of the diffuse gas in recent hydro dynamical cosmological 
simulations that follow the formation of a disk galaxy inside the same halo used in this 
paper 16 . Indirect estimates of the present-day density of the gaseous corona surrounding 
the Milky Way based on OVI,OVII and X-ray absorption measurements are also in the 
same range (~ 10 -4 cm~ 3 at a distance of 50 kpc 42 ). The temperature of the halo is 10 6 
K at 50 kpc, which is consistent with the same observational constraints. Assuming an 
isotropic model, the halo temperature at a given radius r is determined by the cumulative 
mass distribution M(r) of the dark, stellar and gaseous components of the Milky Way 
beyond r and by the density profile Ph(r) of the hot gas 27,43 
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where m p is the proton mass, G and ks are the gravitational and Boltzmann constants. 
In reality, the gaseous halo will not be in hydrostatic equilibrium since it would slowly 
cool and sink to the center while at the same time new gas flows within the virialized 
region from outside. However, cosmological hydrodynamical simulations 16 indicate that 
the characteristic density of the halo at 30 — 50 kpc from the center does not change 
significantly from z = 2 to z = 0. 

5 Tests of the results 

5.1 Tidal heating and gas removal 

We have performed a number of tests to explore in detail the nature of the transformation of 
the baryonic component. The main goal was to determine the main cause of the expansion 
and roundening of the stellar component. We investigated the relative role of the rapid 
removal of the gas, which comprises most of the baryons initially, and that of the repeated 
tidal shocks in changing the shape of the stellar component. Although the initial dwarf 
model is dark matter dominated by construction, baryons comprise as much as 30% of the 
total mass within the initial exponential disk scale length. Their rapid removal due to ram 
pressure stripping might in principle affect the dark matter potential, possibly lowering its 
density as found in previous studies where baryonic stripping occurs as a result of powerful 
supernovae winds 43,44,45 . 

Figure 4 shows principle axis ratios s = b/a and q = c/a (a > b > c) of stars, calculated 
from the eigenvalues of a modified inertia tensor 46 : 1^ = J2a x 'i x 'j/ r ai where xf is the 
i coordinate of the ath particle, r 2 = (y") 2 + (y^/ 3 ) 2 + (2/3 /<?) 2 ! an d yf are coordinates 
with respect to the principle axes. We use an iterative algorithm starting with a spherical 
configuration (a = b = c) and use the results of the previous iteration to define the principle 
axes of the next iteration until the results converge to a fractional difference of 10 -2 . Results 
are shown for three experiments. The comparison between the red and the black lines 
clearly shows that repeated tidal shocks are necessary for the transformation of the stellar 
component to occur; the increase of c/a is amplified by more than a factor of 6 when the 
dwarf is continuously tidally shocked relative to a case in which it just loses most of its 
baryons owing to ram pressure and then evolves without any tidal perturbation (in the 
latter case we removed the primary halo from the simulation after the gas was stripped). 
The green lines in Figure 4 show a test in which about 30% of the initial gas content is 
retained until the end because the UV background is absent and thus the effect of ram 
pressure is reduced (we tested that not including the cosmic UV background at all or 
including it with the low amplitude expected at z < 1 from the Haardt & Madau model 
yields the same result 26 ). The residual gas is funnelled to the center by the bar, where it 
forms a very dense knot 26 . We notice that despite the fact that this is a small fraction of 
the initial gas content it represents a non-negligible contribution to the central potential 
because it all ends up in the inner 0.5 kpc (corresponding to the initial disk scale length). 
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Figure 4. Shape evolution of the stellar component. The upper panel shows the intermediate- 
to-major axis ratio, b/ct, and the bottom panel presents the minor-to- major axis ratio, c/a. Solid 
lines correspond to the self-consistent simulation in which tides, ram pressure stripping, and UV 
background were included throughout the orbital evolution of the dwarf. Dotted lines correspond 
to the case where the dwarf galaxy was extracted from the simulation just after the gas was 
removed and evolved in isolation. Dashed lines show results of a test in which approximately 30% 
of the initial gas mass remains in the galaxy because the absence of the UV background reduces 
significantly the effectiveness of ram pressure stripping. Thick and thin lines correspond to 10 
and 5 Gyr of orbital evolution, respectively (for the extracted dwarf we only show the result after 
10 Gyr). The shape of the stellar component changes gradually over time and the presence of UV 
heating is essential for the shape transformation of the dwarf. 
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In this case the dwarf maintains a bar-like, prolate shape even after 10 Gyr, hence the 
lower c/ a. After 5 Gyr the central density is 8 times higher compared to the simulation in 
which complete gas stripping occurs. As a result £ s hock > ^orb within 0.5 kpc, where t s h oc k 
is .Rperi/Vjperi, -Rperi and V per i being, respectively, the pericentric distance and the velocity 
at pericenter, and t or b is the orbital time at 0.5 kpc from the center of the dwarf; the 
response of the system to the tidal forcing is thus adiabatic instead of impulsive, hence 
further morphological changes are inhibited. This demonstrates that while gas removal 
does not drive the transformation it represents a crucial step for the evolution of the stellar 
component because the effect of the tidal shocks is considerably weakened when a dense 
central gas component is retained. 

Figure 4 also shows that c/a does not grow beyond 0.2 in isolation, even after a few Gyr 
of evolution. The average (apparent) ellipticity measured for late-type dwarfs in the Local 
Group is larger, ~ 0.58, for the 8 dwarfs having a range of rotational velocities comparable 
with those of our initial dwarf galaxy model (Vp ea k ~ 30 — 40 km/s) 47 . Similar values 
are found for dwarf irregular galaxies in clusters 48 . However, in comparing with observed 
values of ellipticity one must take into account inclination effects, since the c/a discussed 
for our simulated dwarf is the intrinsic value, i.e. the value that would be measured by an 
observer viewing the dwarf exactly edge-on. By viewing the dwarf at inclinations in the 
range 30-60 degrees, which statistically are a much more representative situation, the mean 
apparent c/a varies in the range ~ 0.4 — 0.7 within 3 disk scale lengths (it is 0.55 at 45 
degrees of inclination), hence perfectly consistent with the observed values. 

The central dark matter density decreases by almost a factor of 4 during the evolution 
as a result of the repeated tidal shocks when all the gas is readily removed. Since the dwarf 
is dark matter dominated, it is natural to ask whether the change in the halo potential 
associated to this drop of the halo central density is sufficient to alter the mass distribution 
of the stellar component or if the direct heating of the disk by the tidal shocks is crucial. 
Figure 5 shows thickness of the stellar disk defined as the dispersion of particle positions 
above the midplane (z = 0), as a function of cylindrical radius R. The black line corresponds 
to the initial disk and red line corresponds to the disk thickness after 5 Gyr of evolution in 
isolation. The good agreement between these two curves suggests that the initial disk model 
is in reasonable equilibrium. The green line shows instead the results of an experiment with 
the isolated dwarf galaxy model in which we decreased the mass of the inner dark matter 
halo over 5 Gyr by the same amount as in the simulation in which the primary halo is 
included. The decrease of the central dark matter mass is of order of 40% over 5 Gyr and is 
responsible for an increase in thickness of less than 30%. The stellar component maintains 
an evident disk-like structure. Hence the transformation of the disk is not driven by the 
response to the changing dark matter potential. 

In summary, these experiments indicate that the repeated tidal shocks in the impulsive 
regime are driving the transformation of the shape and mass distribution of the dwarf 
galaxy and that gas removal is a necessary condition for them to be effective. We note that 
in previous works on tidal stirring 49 ' 50 morphological changes were happening primarily as 
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Figure 5. Evolution of the disk thickness as a function of radius, computed relative to the initial 
disk scale height, that characterizes the vertical surface density profile (Rd is the initial disk 
scale length). The results are shown for the run in isolation (red line) and for the run in which 
40% of the halo mass within a sphere containing the disk radius is gradually removed over a 
timescale of 5 Gyr (green line). The initial disk thickness is also shown for comparison (black 
line). 

a result of tidally induced dynamical instabilities rather than impulsive heating. There 
the buckling of the bar was driving the change of shape from a disk into a spheroid. Here 
buckling never occurs because the self-gravity of the stars is too low for global bending 
modes to grow. In fact in those previous works the dwarf models had similar baryonic 
masses but stellar masses nearly an order of magnitude larger (conversely the gas was at 
most contributing 30% of the baryonic mass). 

Although tides and non-axisymmetric instabilities heat the system and remove angular 
momentum 49,50 , some residual rotation is present in the dwarf after ten billion years. Figure 
6 shows the intrinsic rotational velocity as a function of radius; this corresponds to the 
maximum rotational velocity that can be observed, i.e. the velocity measured by observing 
the dwarf perpendicular to its rotation axis. The rotation is negligible near the center but 
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Figure 6. Rotational velocity profile of the stars as a function of the distance from the center 
of the dwarf after 10 Gyr of evolution. The profile is azimuthally averaged and spherical radial 
bins have been used. 

becomes larger than ~ 3 km/s at several hundred parsecs from the center. Interestingly, 
this is comparable to the rotation inferred for Ursa Minor 51,52 , for which a velocity gradient 
around the apparent morphological major axis has been observed out to a distance of a few 
hundred parsecs. How much such rotation reflects the original rotation of the dwarf and 
how much it is the result of dynamical evolution driven, for example, by tidal torques, will 
be subject of a future detailed investigation. 

5.2 Resolution tests on gas density 

Our simulations lack radiative transfer. Self-shielding of the gas could reduce the net heat- 
ing from the cosmic UV background and therefore reduce the effect of ram pressure. Re- 
cently, cosmological simulations with SPH hydrodynamics have been used to study the 
formation of dwarf galaxies at high redshift, before as well as during the reionization 
epoch 53 . These include radiative transfer to follow the propagation of UV photons pro- 
duced by a cosmic ionizing background similar to the one employed here and find that 
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the gas self-shields efficiently only for densities p > 10~ 22 g/cm 3 . The gas has always a 
density at least an order of magnitude lower than the latter threshold in our simulations; 
even during the phase in which the density increases the most, namely after the bar-driven 
gas inflow, we measure p ~ 10~ 23 g/cm 3 . We conclude that self-shielding will be negligible 
given the low density of the gas in the simulated dwarf. However, resolving the highest gas 
densities correctly is not trivial in numerical simulations. We tested that our inference on 
the importance of self-shielding is not dependent on the adopted resolution by running the 
same simulation at three different mass and force resolutions, spanning a factor of 64 in 
number of particles (the highest resolution run employs 1.6 million gas particles and was 
run only for one orbit). The results are shown in Figure 7. As it can be seen, between the 
lowest and the highest resolution case the difference in the maximum gas density is less 
than a factor of 2, and it is less than 40% between the intermediate and highest resolution 
cases. The intermediate resolution case corresponds to the standard resolution adopted in 
our simulations. The figure also shows that a decrease of a factor of 3 in the gravitational 
softening implies an increase in the maximum density of only about 30% at fixed particle 
number. 

6 Critical discussion of the results 

6.1 The various modes of ram pressure stripping 

The gas is first stripped by tides and instantaneous ram pressure as the galaxy approaches 
pericenter. This removes the outer disk, outside about two initial disk scale lengths. For 
such gas, the ram pressure force, p g h p V P 2 (Pgh p is the density of the halo near pericenter 
and V p is the orbital velocity of the galaxy near pericenter, that are, respectively, 10~ 4 
atoms cm -3 and 300 km/s) is stronger than the gravitational restoring force, which is 
mainly provided by the dark halo. More instantaneous ram pressure stripping occurs also 
at second pericenter as the halo central density, and hence the gravitational restoring force, 
has decreased in the meantime due to the tidal shocks, and the galaxy loses its gas down to 
about one initial disk scale length. Quantitative comparisons with analytic predictions 54 for 
instantaneous stripping are difficult because tidal forces continuously affect the magnitude 
of the gravitational restoring force and participate in the stripping, and the gas is subject 
to compressional heating, radiative cooling and an external radiation field. Nevertheless, it 
is important to notice, that dwarf galaxy models almost identical to the one simulated here 
have stripping radii within 30% of those predicted by the Gunn & Gott formula 54 when 
they are evolved in periodic tube flows that simulate a single pericenter passage in absence 
of tidal forces 26 ' 42 . 

However, the gas inside about one disk scale length lies inside the tidal radius at all 
times, hence tidal stripping does not play a role; in addition, for such gas the gravitational 
restoring force is always higher than the instantaneous ram pressure force. Stripping might 
still occur by ablation via Kelvin-Helmoltz instabilities, the so called turbulent stripping, 
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Fi gure 7. The maximum density of the gas for runs with different mass and force resolution. For 
all the runs the density is measured after bar formation has triggered the gas inflow and its value 
relative to that in the run employed in the Letter is shown (p\ = 2 x 1CP 23 g/cm 3 ). Red open 
circles refer to runs in which the gravitational softening is equal to the run discussed in the Letter 
(35 pc), while black filled circles show the results for runs in which the softening was increased 
by a factor of 3. Therefore, while red and black circles represent the results of two sequences of 
runs with varying mass resolution, comparing red and black circles at fixed number of particles 
shows the effect of changing the gravitational softening. 

or by laminar viscous stripping 26,37,55,56,57 . As we have also shown elsewhere 26 , Keivin- 
Helmoltz instabilities should not occur for dwarfs embedded in a massive dark halo because 
of the stabilizing effect of the gravity of the halo. This is the case even here. Following 
existing analytical predictions 55 we estimate that the critical value of the gravity (in sim- 
ulation units) to stabilize our gas disk against Kelvin-Helmoltz is g cr ~ 0.1 initially, where 
g cr = 2nV 2 /DR g (V ~ 200 km/s is the average orbital velocity of the galaxy, R g = 0.5 
kpc is the scale length of the gas component, and D = p g /p g h is the ratio between the 
density of the gaseous disk at R ~ R g and the time averaged density of gaseous halo along 
the orbit of the galaxy). The self-gravity of the dwarf galaxy within R g , which is mainly 
provided by the dark halo, is instead g ~ 2.5 > g cr , where g = GM ha [ j R g . We note that 
g > g cr holds down to about 1/5 of the disk scale length, and also for R> R g . Hence the 
gaseous disk of the dwarf is initially stable to Kelvin-Helmoltz instabilities. 
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When g > g cr the gas can still be stripped by laminar viscous stripping. The mass 
loss rate in this case is given by M ~ 12/2.87rR g 2 V(X max /R g )(ci/V) 55,5 , where A max is the 
longest wavelength not stabilized by gravity. Application of this formula to our dwarf galaxy 
in motion through the hot galactic halo yields a characteristic timescale for laminar viscous 
stripping which is longer than the Hubble time, mostly because A max < R g . However, close 
to the second pericenter passage, the situation changes (see Figure 8). The density of the 
halo within 0.5 kpc decreases by a factor of 4 as a result of the tidal shocks, so that g drops 
to ~ 0.6 at one disk scale length, now very close to g cr (the latter has also decreased in the 
meantime, g cr ~ 0.3, since the mean density of the gas has dropped). The gas component 
of the dwarf is now only marginally stable to Kelvin-Helmoltz instabilities, but using the 
formula to calculate the mass loss rate due to Kelvin-Helmoltz 55 for R g = 0.5 kpc one 
derives a characteristic stripping timescale still quite long, i.e. comparable to the Hubble 
time. Indeed, in the simulations we do not see Kelvin-Helmoltz instabilities developing 
before the gas is removed. Instead, assuming that the longest unstable wavelength is 
Amax ~ 0.5 kpc one obtains a laminar viscous timescale of ~ 2 Gyr. The gas is indeed 
completely stripped ~ 1 Gyr later (it is possible that stripping is slightly enhanced by 
artificial viscosity of SPH 26 ). 

When the UV background is absent g » g cr even at the second pericenter passage 
because (1) the density of the halo within 0.5 kpc decreases by only 50% as the cold, 
concentrated gas reduces tidal heating (see section 5.1 above), and (2) g cr is smaller because 
the gas is colder and denser relative to the case where the UV radiation background is 
included, thereby increasing the value of D (see above) by a factor of 3. As a result, the 
system is very stable to disruption by Kelvin-Helmoltz instabilities. Likewise, the longest 
unstable wavelength is much smaller than the radius of the residual gas, which implies a 
laminar viscous timescale close to the Hubble time, similar to what we have obtained for the 
initial conditions. This is consistent with our result that the gas is not completely stripped 
when the UV background is absent. In conclusion, despite the strongly time-dependent 
behavior of all the variables involved in the problem, the results of the simulations are 
consistent with the expected characteristic timescales and relative importance of the various 
ram pressure stripping modes. 

6.2 Adding complexity to the scenario 

The "proximity" effect due to the local ionizing background from the Milky Way at z > 1, 
not included here, could have been even more important than that of the uniform cosmic 
UV background. Galaxy formation simulations that include star formation and the effect 
of supernova feedback show that star formation rates comparable to those of a luminous 
infrared galaxy (LIRG), > 10M o /yr, are attained around z ~ 2 for a Milky Way-type 
galaxy 16 . We have run a test with the addition of the local ionizing radiation produced by 
the starburst, found complete ram pressure stripping already on the first orbit and, as a 
consequence, an even faster transformation of the stellar component. 
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Figure 8. Color-coded projected density map of the gas component. The box is 3 kpc on a 
side and shows the gas distribution view from above the orbital plane. The time is close to the 
second pericenter passage, after about 3 Gyr of evolution. The box corresponds to a zoom in of 
the middle right panel in Figure 1 of the Letter, the density limits are also similar, from 1CP 31 
g/cm 3 (black) to 10 -23 g/cm 3 (yellow). 

The processes discussed in this paper assume that dSphs are bound satellites of a 
primary galaxy. Isolated dSphs such as Tucana are found at distances of order twice the 
expected virial radius of the Milky Way or M31 47,50 . Either they are not bound to either 
of the primary galaxies or are extreme satellites on nearly radial orbits (in cosmological 
simulations a few satellites whose orbital apocenter greatly exceeds the virial radius are 
typically found 18,47 ). In either case they cannot be objects that fell early into the halo 
of a primary galaxy (if bound, their orbital time would be too long to be consistent with 
an object infalling at z > 1). Instead, their mass-to-light ratios are likely dictated by 
their halo mass rather than by the influence of the environment. Either their halos have a 
relatively small mass, V^eak < 20 km/s, so that reionization alone could have been sufficient 
at suppressing their baryonic content and star formation by inhibiting the collapse of gas 
and/or photoevaporating already accreted gas 53 , or their halos have masses comparable 
to those of Draco but their mass-to-light ratios are much lower, comparable to those of 
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the other dSphs, because baryonic sweeping was not effective. Transformation from disk- 
like progenitors into dSphs could still be the origin of the diffuse spheroidal structure of 
Tucana since tidal stirring can occur even outside the virial radius of the primary galaxies 
in cold dark matter models, for example between members of subgroups infalling into the 
primaries 58 . Improving our knowledge on the structure and properties of isolated dSphs 
can thus yield important constraints on formation scenarios such as that proposed in this 
Letter. Obtaining internal velocities to measure their halo mass and determine their mass- 
to-light ratio as well as achieving accurate upper limits on their gas content should be 
priorities among the new observational programs of nearby dwarfs. 

Finally, we have only considered NFW profiles for the halos of the dwarfs. Recent 
work 45 has shown that a cuspy dark matter halo fitted by an NFW or Moore profile might 
develop a core as a result of a strong mass outflow such as that caused by supernovae winds 
after a burst of star formation. Our scenario for baryonic sweeping would be even stronger 
if dSphs had cored halos. Indeed, as we have seen in the previous section, the various 
modes of ram pressure stripping are all enhanced if self-gravity is lower, as it would be the 
case for a cored halo as opposed to a cuspy halo. Stripping will be much faster especially 
if the size of the core is comparable to the disk scale length since, as we have explained in 
section 6.1, the inner region is that in which instantaneous stripping is not effective and 
other slower stripping regimes such as viscous or turbulent stripping come into play. Tidal 
shocks will also be enhanced because the response of the system will be more impulsive, an 
effect which is already seen when comparing subhalos with NFW profiles having different 
concentration parameters 59 . On the other hand, dwarf galaxies having the low gas surface 
densities and high Toomre parameter of our initial model would hardly undergo a strong 
burst but rather have a fairly low star formation rate 35 , in which case modifications of the 
halo profile are unlikely. 
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